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Abstract 

The gravitational instability of Yang-Mills cosmologies is numeri- 
cally studied with the hamiltonian formulation of the spherically sym- 
metric Einstein- Yang-Mills equations with SU(2) gauge group. On the 
short term, the expansion dilutes the energy densities of the Yang- 
Mills fluctuations due to their conformal invariance. In this early 
regime, the gauge potentials appear oscillating quietly in an interac- 
tion potential quite similar to the one of the homogeneous case. How- 
ever, on the long term, the expansion finally becomes significantly 
inhomogeneous and no more mimics a conformal transformation of 
the metric. Thereafter, the Yang-Mills fluctuations enter a complex 
non-linear regime, accompanied by diffusion, while their associated 
energy contrasts grow. 



PACS numbers: 04.25.Dm, 98.80.Jk 



1 Introduction 



The spherically symmetric Einstein- Yang-Mills (EYM) equations with 577(2) 
gauge group has been extensively studied in the past decades, from cosmolog- 
ical solutions (cf. [H El El El E] ) to static configurations (cf. jHlIZI)- The ho- 
mogeneous and anisotropic cosmological models with Yang- Mills (YM) fields 
were considered in [HJ El HD] and exhibit chaotic features. For a detailed 
review of the litterature prior to 1999, see [Zj and references therein. More 
recently, some authors have investigated non-abelian Born-Infeld cosmology 
jTTj in which the YM-Born-Infeld lagrangian is no more scale invariant and 
therefore yields to different behaviour from classical YM fields, especially 
near the singularity. 

However, most of the work done so far provides (semi-)analytical and/or 
numerical solutions in the case of only one independent variable: time (cos- 
mological solutions) or radius (gravitating non abelian solitons and black 
holes). Nevertheless, some authors [5] have proposed time-dependent in- 
homogeneous analytical solutions of EYM theory but at the cost of strong 
restrictions on the form of the metric and the gauge potentials. As well, per- 
turbation methods have been used several times to study the instability of 
self-gravitating non-abelian solutions [?.. Here, we use the hamiltonian for- 
mulation of spherically symmetric EYM equations with SU (2) gauge group 
to build a numerical method that allows us to find numerical time-dependent 
inhomogeneous solutions under a few general assumptions. In this paper, we 
apply it to the question of gravitational instability of YM cosmologies, by 
studying collapsing shells, which was done for Minkowski space in [TJ|. 

In the case of homogeneous and isotropic spacetimes, the conformal invari- 
ance of YM fields allows to solve separately Einstein and YM equations as 
the resulting solution is a radiation dominated universe. The general ansatz 
for open, flat and closed cosmologies filled with su{2)— valued YM fields was 
proposed in [I], for which the solution in is a special case. In this paper, we 
will investigate departures from this ansatz, consequent to some initial per- 
turbations, in order to characterize their gravitational instability. Although, 
our approach differs from usual perturbation theory in the sense that we 
do not linearize the equations, we consider here time- dependent inhomoge- 
neous numerical solutions around the general ansatz. One will argue that 
the conformal invariance of YM fields implies that any primeval excitation 
of these fields will be diluted as the universe expands. We will see that this 
is true only in the first times of the evolution of the fluctuations, when the 
expansion is not too inhomogeneous and looks like a conformal transforma- 
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tion of the metric. Thereafter, interesting phenomena such as diffusion and 
self-interaction arise which results in growing instabilities. 

In section 2, we remind the reader about the hamiltonian formulation of 
the spherically symmetric EYM system with SU(2) gauge group, firstly in- 
troduced in then we recall the Gal'tsov-Volkov ansatz that we will use 
as a background solution. Thereafter, the numerical method used to inte- 
grate the hamiltonian formulation is introduced in section 3, a test on the 
homogeneous ansatz being provided. The section 4 is devoted to a quali- 
tative discussion of the gravitational instability from the numerical results 
presented in the figures at the end of the paper. Finally, we conclude by 
some comments on the application of our results to cosmology as well as 
some perspectives. 

2 Theoretical Framework 
2.1 Hamiltonian Formulation 

The hamiltonian approach of the gravitational field equations is the well- 
known Arnowitt-Deser-Misner formalism which can be completed by 
a hamiltonian treatment of YM equations describing a gauge interaction 
|T3|ll5j. The spherically symmetric inhomogeneous vacuum space-times were 
studied in [TT)|IT7] while the spherically symmetric EYM equations with gauge 
group SU(2) were first discussed in |13j . We briefly recall here the hamilto- 
nian equations of the EYM system and fix the gauges in which we will work 
for the rest of this paper. 
The coupled EYM action can be written as 

S = JyTl - *x (1) 

where k = 8irG (c = 1), g is the determinant of the 4-metric, 1Z is the scalar 
curvature, F* v are the components of the YM field strength tensor 1 . For 
a particular gauge group G locally specified by its Lie algebra [T a , Tb] = 
if abcTc, the field strength tensor is given through its components by F° = 
d^A^ — d v A c + fabcA^A^ where the A*'s are the components of the one-form 
connection A = A^dx 11 = A*T a dx^- Here, we will focus on the gauge group 
SU(2) with generators T a = |r a where the r a are the usual Pauli matrices 
and the structure constants are given by the completely antisymmetric tensor 

1 Gauge indices will be indicated by bold latin letters or bold numbers while latin indices 
are only spatial and greek indices are for spacetime. 
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/abc = e abc- The normalization is Tr(T a T h ) = |<5 a t>- 

Assuming spherical symmetry of the spatial slices of the universe, we write 
the metric in a 3+1 decomposition as follows : 

ds 2 = (-N 2 + N X N X ) dt 2 + 2N x dxdt + e 2 ^d X 2 + e 2X dtt 2 , (2) 

where dVL 2 = d6 2 + sin 2 9 dip 2 is the solid angle element, N and N x are the 
lapse and shift function respectively 2 , all these components N, N x , /i, A 
being functions of the coordinates x an d t. 

If the gauge fields exhibit spherical symmetry up to a gauge (see [T3J H3 
Ej), they can be written thanks to the so-called Witten ansatz JH] (with the 
same conventions as in 

A^dx* 1 = aT 3 dt + bT 3 dx + (cTi + dT 2 ) dd + 

(cot 6 T 3 + cT 2 — dTi) sin 9d(p, (3) 

where a, b, c, d are all functions of both coordinates \ an d t. In Eq.@, 
only three components are of physical relevance since the potential d can be 
gauged away by a transformation that preserves the SO (3) symmetry (see 
[T| ITT] and references therein for a detailed discussion). 
If we now let the spatial components of the metric and the connection Af 
(i = 1, 2, 3) be canonical variables and if we consider the lapse and shift 
functions N, N x as well as the electric components Aq of the YM connection 
as being Lagrange multipliers, it is possible to define some momenta conju- 
gated to these variables and then to rewrite the coupled EYM action Eq.((TJ) 
under constrained hamiltonian form (cf. ^SIEIE] and references therein): 



S 




{<K ij d g tJ + doA, a - NH ± - N*Hi - A*g a ] dt d 3 



(4) 



where = \/^g (g % ^K\ — K 1 ^ (K % i being the extrinsic curvature tensor), 

/(3) 

t^a = n 9 (d^Foi ~ N k g tj F* i ) . In Eq.(JIJ), the generators of normal and tan- 
gential deformations (H± and Hi, respectively) and those of internal gauge 
transformations Q a are given by 3 



H 



-fj 



R + 9 1 1 



TCi 



/ 7 a 7 
-9ij iy A 71 A a 



B i a i3 i a ) « 



k a 



—IT 



A aA 



i + f. 



ab^A c 



A) 



2 The angular components of the shift vector Ng - 
due to spherical symmetry. 

3 The gauge coupling constant has been set to one. 




» 

g Q2 and N v 



(5) 

(6) 
(7) 

<7o3 vanish exactly 



4 



with R the scalar curvature of the 3-metric and tt = 7r£, | denoting the 
covariant derivative w.r.t. the 3-metric and the vector B l a reading 

K = \& k (4a,i - A ja>k - f cha A c 3 A^ • (8) 

Let us now focus on the spherically symmetric hamiltonian formalism. 

For the gravitational part of the EYM system, the canonical variables are 

chosen as 4 gij = diag (e 2/ \ e 2A , e 2A ) while their conjugate momenta are given 

by 7Ty = diag (e- 2 ^,e~ 2X ^,e~ 2X ^y Variation of Eq.® w.r.t. the 

lapse and shift functions N and N x will provide the hamiltonian and super- 
momentum constraints while the variation w.r.t. the canonical variables fi, 
A and their conjugate momenta 7r M and it\ will lead to the Hamilton equa- 
tions. As we are interested in perturbations of homogeneous solutions and 
as the YM lagrangian in Eq.(£Q) is conformally invariant, we will work in the 
conformal gauge for the gravitational field: N(x, t) = R(t) where R(t) is 
the scale factor given by the Friedmann equation and N x (x,t) = 0. After 
the variation and the choice of gauge, we now have for the dynamics of the 
gravitational field the following constraints: 

Hi = = -/j'tt^ - A'tta + tt^ - 167riVe 2A+/x 7? = (9) 

H = = ^- 2A {8e 4A (-VA , + 6A' 2 + 4A")- 16e 2 " +2A + vrj - 27^} 
+ 167re 2A+ ^ T ° = 0, (10) 



and the following Hamilton equations, 

/'< 

A 

^A 



(11) 

AV A )+7T M (7T M -27r A )} 
.A' 2 -A") + ^(^-2tta)} 



= \e-»- 2X N(n, - tta) 

-16tt N e 2A+/i Tl 
-32tt N e 2X+fl Tl 



where are the components of the YM stress-energy tensor (see further), 
and where a prime and a dot denote derivations w.r.t. position x an d time 

4 We work in the following non coordinate basis : (d x ,dg, ^-§d,p) 
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t respectively. 

On the YM side, we have, according to Eq.©, only one electrical component 
Aq = a and three couples of variables and momenta (Af, vr^' 3 ) = (b, 7Tb), 
(Al nf) = (Al nf) = (c, vr c ), (Aj, nf) = -{A*, nf) = (d, n d )- 
Varying the YM action w.r.t. these variables and fixing the gauge d(x, t) = 

so that iTd = j^ ac leads to the following constraint 

g 3 = = ix' b + 2^ac 2 = 0, (12) 

and to the following Hamilton equations 

b = a' + Ne^- 2X n b 

c = 7Ve _/ V c 

7f 6 = -2Ne~^bc 2 (13) 

7f c = ^a 2 c + (Ne-^)' - Ne-% 2 c - Ne^ 2X c(c 2 - 1) 



N 2 e" 2 " ( b ' + 2 b ^-^b]-2Ne-^-^e 
\ c J c 4 



/J.-2X 



71^ - 7T X )a 



N 
+ N a - 



With these variables, we can now write the components of the YM stress- 
energy tensor = -2F* a F u aa + \g^ u F^F a ^ 



it 2 r' 2 N 2 + p 2 ^n 2 r 2 4- N 2 b 2 r 2 

2 ~ e 2 c N 2 

+e" 4A ^£ (14) 

T i ,-4A ^ _ - 2 A- 2/V 2 _ c 2 iV 2 + eVc 2 + Wc 2 A _ 2fl 

2 1 - e 2 c N 2 

+e" 4A ^£ (15) 

-4A 

Tl = -—(*l + {c 2 -l) 2 ) (16) 

-2A 

T i° = 2 ^ (^e^vrec' + abc 2 ) = -N 2 e~ 2 ^, (17) 

in order to achieve the coupling between gauge and gravitation sectors. 
Eqs. (JHl fTTj) are the hamiltonian formulation of the well-known spherically 
symmetric EYM equations with gauge group SU(2). In the next paragraph, 
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we recall the solution for those fields in a Friedmann-Lemaitre universe be- 
fore describing the method used to study gravitational instability of such 
solutions. 



2.2 Yang-Mills Cosmologies 

Non—abelian YM fields can fill a non trivial, radiation dominated, Friedmann- 
Lemaitre universe due to their increased number of degrees of freedom, and 
the couplings between them, compared to the abelian Maxwell theory. In 
particular, the case of gauge group SU{2) has been extensively studied, so 
we refer the reader to [TJ |3J 12 HI HH Hj and others references therein for a 
review. 

With the definitions of previous paragraph, we can write the homogeneous so- 
lutions we will use for the background in which the fluctuations will evolve. 
First, the geometry is given by (in the conformal gauge N = R(t), upper 
script B stands for background): 

H B = logR 

X B = lo gj R + log£ 

vrj = —ARRTj 2 (18) 
Trf = -8RRY, 2 

where £ = sin x, X-, srn h X f° r closed, flat and open cosmologies respectively. 
The conformal invariance and the consequent vanishing trace of the YM 
stress-energy tensor result in a radiation-dominated universe for which the 
Friedmann equation Eq.([TU| reads (using Eqs. (fT%|) ): 

H 2 =\^ \ =^"*. (19) 

with £ some constant to be precised in one moment. This allows us to write 
the following solutions for the scale factor R(t) (see also [211 )• 

R(t) = J^Ct(k = o) 

R(t) = J^C cos(t) (k = +1) (20) 

R(t) = W~C sinh(t) (k = -1) 
V *^ 

for flat, elliptic and hyperbolic universes, respectively. The constant C = pR 4 
is a first integral of the conservation equations V M T^ = 0. 
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The gauge fields and their conjugate momenta are given by the following 
ansatz [TJ [TT] : 



c B ( X ,t) = y/l + E 2 (a 2 - k) 



1 + E 2 (a 2 - fc) 

(^ 2 - k) 
1 + E 2 (a 2 - k) 



b B (x,t) = aE 2 ,fr 2 , - (21) 



R . aaE 2 
vr c B (x,t) 



y/l + E 2 (a 2 - k) 

where k = +1, 0, —1 for closed, flat and open cosmologies respectively, 
and a is a function of time. By substituting the ansatz for the gauge fields 
Eqs.(j2U) Eqs. ([12M13p . the YM equations now reduce to 

a + 2a (a 2 - k) = 0, (22) 

or, expressed as a first integral 

d 2 +(a 2 -A;) 2 = £ 2 , (23) 

where S 2 = |C = |pi? 4 . The conformal invariance of YM fields is now man- 
ifest as the scale factor R(t) does not appear in Eqs. (j22p23j) . Therefore, in 
the conformal gauge, the time part of the gauge fields a(t) oscillates in the 
potential well (a 2 — k) 2 . In the synchronous gauge N = 1, it can be shown 
that the frequency of the oscillation decreases with expansion. 
In the rest of the paper, we will focus on the case of flat space as the ob- 
served features of gravitational instability can be deduced from the general 
properties of the homogeneous solution presented in this paragraph and can 
be therefore extrapolated to curved cosmologies. 



3 Overview of the Method 

In order to integrate the hamiltonian EYM equations, we will use a two- 
step numerical procedure often referred as free evolution method. First, we 
solve the three constraints Eqs.([§lfTUj) and Eq.([12|) - this consists of the ini- 
tial value problem (IVP) - then we propagate the initial data of the fields 
imposed or found by the IVP with the Hamilton equations Eqs. pillT^j) and 
check the accuracy and stability of evolution by evaluating the error made on 
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the constraints. The growing violation of the constraints in such procedures 
and the development of formalisms reducing it are very hot topics in numer- 
ical relativity (see for example [20] , and j2l] for an introduction to numerical 
relativity). Here, we choose the canonical approach to the EYM system (orig- 
inal ADM formalism ^3] for gravitation and classical treatment for YM [T3*] ) 
which we will rewrite in section 3.2 in order to improve stability and accuracy. 

But before going further, let us give some comments on the particular choice 
of gauges we made for both the gravitational and YM fields. For the first 
one, we chose a slightly modificated version of the famous geodesic slic- 
ing: N(x,t) = 1, Ni(x,t) = 0, which we called the conformal gauge: 
N(x,t) — R{t) and N t (x, t) = 0. Even though the geodesic condition is 
well-known to produce coordinate singularities (cf. |22| ) . it is the most ef- 
ficient for weak gravitational fields such as instabilities in a homogeneous 
background. The reason why we chose a conformal gauge with iV = R(t) 
rather than a synchronous one N = 1 is because the YM fields oscillate at 
a fixed frequency in the first while their period is stretched with time in the 
second. Therefore, using a synchronous gauge would have increased the time 
over which the fluctuations around the homogeneous solution significantly 
change. And, as we used an explicit scheme, it was crucial to reduce this 
time in order to avoid numerical discrepancies. 

Furthermore, our choice of a gauge dependent formalism (d was set to 0) 
was prefered for two reasons. The first is that it will sweep away one degree 
of freedom in the hamiltonian formalism, which simplifies the equations and 
therefore lowers truncation errors as well as avoiding those related to the cou- 
ple (d, 7Td) of unknowns. The second is because it simplifies the algorithm. 
Indeed, with a gauge-invariant formalism, the electric potential A$ = a which 
represents the freedom of making arbitrary gauge transformations in the YM 
fields during the evolution of the system has to be determined by an addi- 
tionnal gauge condition. By fixing the gauge, we are left with a Hamilton 
equation for this component. However, by reducing in this way the number 
of degrees of freedom in the analysis, we cancel an additional scalar field in 
the theory, in the context of the two-dimensional abelian Higgs mechanism 
that consitutes the spherically symmetric EYM system with gauge group 
SU{2) (cf. [IBJIIj)- But now, we have learned from a more deep analysis of 
our results that these scalar fields and the interactions related to it bring a 
good part of the energy of the fluctuations. Part of our future work will be 
to examine how energy spread out into the different modes of the fully gauge 
invariant YM equations and how it affects the results presented here. 
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3.1 Initial Value Problem (IVP) 



In this problem, we have to determine the initial distribution of nine fields 
(/io, Ao, vr Mi0 ) tta,o, a o, bo, Co, tt^o, ^c,o) 5 that verifies the three constraints 
(two for the gravitational field and one from the gauge sector) of the EYM 
hamiltonian system. Here, we will only write the formulas for a perturbation 
of a flat Friedmann-Lemartre universe as they can easily be generalized to 
closed and hyperbolic cosmologies. In order to do this, we fix six fields as 
follows: 

• We start working some time after the Big Bang to avoid primeval sin- 
gularities : .R(O) = 1. Please note that the formula given in section 2.2 
for the scale factors should be translated accordingly. 

• We set Xq = lnx so that the initial circonferential distances are taken 
as references for the coordinate x- 

• We perturb the following fields around their homogeneous values : 



where is the initial background density, taken as a parameter. The 
functions (x) are also provided and will be taken so as to have con- 
venient boundary conditions for the fields. Typically, we will study the 
gravitational instability of shells in a Friedmann-Lemaitre background, 
i.e. we will impose that the e^'s vanish at the boundaries of the x co- 
ordinate interval. For example, we can take for the e^'s some gaussian 
functions. 

From there, we solve the super-momentum constraint Eq.Q w.r.t. to ttx,o 
and the YM constraint Eq. (jl2J) w.r.t. /i and we substitute the results into 
the hamiltonian constraint to get the following differential equation for iTb^. 




% +e a (x) 
bo+e b ( X ) 
Co + e c ( X ) 

Tifo + £ c(x) 



(24) 



A4>" + B4>' + C(p' 3 - 3 = 



(25) 



5 Lower script zero means at the time t = 0, for example /io = £t(x, 0)- 
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where we put <fi = 7r bi0 (x) = 7T&(x, 0) and where the coefficients are 
A = 8a 2 4x 3 



B = -A X 2 [2a' a c 4 x + 8irc' 2 a 2 4 + Ac' a 2 c 3 oX + 8na 2 b 2 c 6 + agcgj (26) 

C = 2\' \ x 5 - 8na 2 c 2 x 2 + 167rx 3 ao&o^A - 4ttc^ + 8ttc 2 + 3A^ 4 - 4tt 
+X 2 + 167rx 3 c' AoCo - 8tcx 2 c 2 - 



The fields with lower script mentionned in Eq. (j26|) are those given by hy- 
pothesis in Eq.(jHJ). From the solution (f) = 7r b of Eq. (j2*5|) . one can complete 
the set of initial data with the following relations 

fi = In 



?n 2 r 2 



7T M)0 = -4A e^V 
7T Ai0 = -32ne^c' c x-^ x 3 e' h -32ne^a b c 2 x + 2TT^ (27) 

vr c ,o = e Mo c - 

We used a finite difference scheme to solve Eq. (I25j) as follows: defining \% = 
Xmin + Ax(i — 1) (i — 1, ■ • • , iV x + 1), 0j = <fi(xi) an d analogous definition 
for other fields, Eq. ()25|) becomes a non-linear system of N x — 2 algebraic 
equations for determining the </>j's: 

\ + i - 24>i + 4>i-i , „ - 4>i-i , ~ f 4>i+i - 0i-i x " 



Ax 2 +ft 2A X +t 'V 2A X - 



with the boundary conditions 



0i = flfo(Xi) = vqXmin ( 29 ) 

B / \ 2 



where cr = y — o~q, cr becoming another parameter setting the initial 

position in phase space for the time part of the gauge fields filling the flat 
background. We choose a Newton- Raphson method (cf. j2Sj) to solve Eq. (|2"£J) 
with the initial guess : 0j = n^Xi) = & x 2 and stop the relaxation when 
the error on the equation is saturated. 
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3.2 Propagation of Hamilton Equations 



In this paragraph, we first rewrite the EYM constraints and Hamilton equa- 
tions in order to render the numerical integration more accurate and stable. 
Indeed, a direct integration of Eqs.(JHJ) and Eqs.((THJ) will have to recover 
first the homogeneous solution Eqs. (|18H2~T|) before focusing to its perturba- 
tion. After having first directly integrated Eqs. (fTT|) and Eqs.lfT^j). we found it 
more interesting to rewrite the equations so that they would be more adapted 
to the study of the gravitational instability and less dependent on the homo- 
geneous solution. The change of variables explained below results in a faster 
and more accurate method. 



Let us again assume a flat Friedmann-Lemaitre background and proceed 
to the following change of variables: 

e» = m( X ,t)R(t) 
e 2A = l( X ,t) R 2 (t) X 2 

7T M = -4:R(t) X 2 7T m (x,t) 

tt a = -8 R(t) X 2 Mx,t) 

2 



b(xt) = Pix ' t)x 

[X,) 1 + X 2 7(x,t) 

c(x,t) = V^+xHxJ) 

Mx,t) = np(x,t)x 2 



7Tc(X»*) 



v/l + X 2 7(x,t) 



so that the homogeneous solutions (m B , l B , tc^, n B , a B , j3 B , 7 B , n B , n B ) 

are now a set of constants w.r.t. x '■ (1> 1) R\ Ri o" 3 ; ■, 111 
terms of these new variables, the hamiltonian formulation of EYM system 
Eqs.([niEJ) now become for the constraints 

ml I'y 4- 21 

H t = = — vr mX + -S ^ - 27r - - X< n - 4vr R? m I X 7? = 

ml 

m'R(l' X + 2l) R{l'x + 2lf R(l" X 2 + 4xl' + 2l) rnR 
n = = = : — r - ^ r 



m 2 x 4mlx 2 mx 2 X 2 

- (vr m - 47T;) 

= -n'pX + 2tv/3 — 2ma 
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and for the Hamilton equations 

R 7T m - 

m = ——m — 

R Rl 

1 - -4> +2 ^k < 32 > 

R mR (l' X + 2lfR tt 
* m = -R 7Vm ~W + 8mZx 2 '^dR i7rm ~ A7Tl) 
+4tt R 3 m I T\ 

R m'R(l' X + 2l) R ( 2 (l'x + 2iy 

~2^1R (7Fm " +47lR3ml T ' 
with the following definitions for the components of the stress-energy tensor 

, 2 



1 



7T, 



?x 2 + (7+fy) +p 2 x 2 



2/ 2 i? 4 v " ' ' m 2 I i? 4 (1 + x 2 7) 



a 2 



+ /i? 4 (l + X 2 7) (33) 

1 / i?v (i + x 2 7) v 7+ 2 7 y p y 

(Tj 1 is the same expression as Tq except for the sign of the last two terms). 
Finally, the YM equations can be written 

a' a mnp a (27 + Xl') m 

x x 2 ix 2 (1 + x 2 i) 1 71(51 



(34) 

W - I 3 " 2 ) X 2 ma 2 m! (2 7 + xi) 2 7 + 4^7' + l"x 2 

7T — — |— — — |— 

m (1 + x 2 7) (1 + X 2 7) 2m 2 x 2mx 2 



$ = 


2,9vr 7X 2 


m(l + x 2 7) 


7 = 


2^ 




m 




-2I 


7T> = 


m 



(27 + X7') 2 m m 



7 - — 7 



2 



4m (1 + x 2 7) ^ X 2 I 

a = - 2(3 + (3 ' X + ^Px + a( 7Tm - 27ll + - 
m 2 m 3 \ mlR R 
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In order to integrate Eqs.(j22HSH), we use d a second order difference scheme 
in both x an d t (leapfrog-like) and, for the first step only, a forward time 
centered space scheme (FTCS, cf. j2Hj)- This explicit scheme is easy to im- 
plement but yields a very strong Courant condition in this case. Anyway, for 
the study of gravitational instability during a few oscillations of a in Eq. (|22j) . 
it will be sufficient. 

Another crucial point is the question of boundary conditions, already men- 
tionned in the previous paragraph. As we focus on shells and all our initial 
perturbations have been chosen to vanish rapidly on the edges of the consid- 
ered interval in x, we impose that the boundary conditions follow the homo- 
geneous evolution. For example, a field f(x,t) = fij will have as boundary 
conditions 

h,j = /£ (35) 

fN x +l,j = /jV x +l,j 

where the upper script B stands for the background solution Eqs. (jl8M2T|) . 
3.3 Test on the homogeneous solution 

In this paragraph, we briefly discuss the precision of the algorithm presented 
before. Figures Q to El show the evolution of the mean quadratic error 

JV x i=l 

on the various fields used in the computation. The general divergent be- 
haviour is due to the simple explicit scheme we have used. Therefore, the 
scheme is unstable, partially due to the violation of the boundary conditions 
that introduce truncation errors on the edges of the grid (see next section). 
We found a very strong Courant condition, for example Ax ~ 10 _1 ; At ~ 
10 -6 to have an acceptable precision on the time scale of an oscillation pe- 
riod of the function a B (measuring the period of the gauge fields). For the 
gauge variables a, (3 and 7, we see in Figure [T] that E a > Ep > E 7 , while 
for the metric variables E m > E\ (Figure EJ). For the constraints (Figure EJ), 
the boundary effects dominate the quadratic error on the hamiltonian and 
super-momentum constraints 7i and Tii though the errors are much lower 
away from the edges of the x~ interval (see also Figure ITo^l. These boundary 
effects are present in all the fields in the computation (see also other 3D 
plots). In Figure El (Hi), we see that the error on the gauge constraint Q 3 is 
much lower than those of Einstein's part. 
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4 Discussion 



In this section, we use the numerical results obtained from the method intro- 
duced before in order to present some interesting features of the gravitational 
instability of YM cosmologies. For the sake of simplicity, we used for the ini- 
tial pertubations a gaussian function of the form: 



where e i: xf and w are parameters indicating respectively the amplitude, 
position and width of the fluctuation. The last is related to the coherence 
length Xp of the fluctuation, the length over which the perturbation changes 
significantly. The figures were computed with only one field initially per- 
turbed amongst the set (a , b , Co, Co, ^o). Perturbing several fields will 
change the amplitude of density and pressures fluctuations on the initial 
space-like slice specified by the IVP and will affect the developments of the 
resulting instabilities in the long term. Indeed, as we will see further, the first 
stages of the evolution are characterized by a conformally-invariant regime 
maintained by the too less inhomogeneous expansion. The instabilities arise 
only after some time, depending on the initial conditions, when the fields 
reach a sufficient amount of inhomogeneity. 

Figures El HI M EH an d [T3] show a typical evolution for the metric, the 
gauge potentials and the YM stress-energy tensor components as well as the 
corresponding values of the three constraints under an initial perturbation 
of the gauge field c. It is useful to notice that only 50 iterations have been 
plotted in each 3D plot. 

For the gravitational sector, the perturbations of the spatial coefficients of 
the metric grow up as a power of conformal time (see Figure EJ). This is 
typical of a radiation dominated universe and has been studied for a long 
time by several authors such as Lemaitre, Tolman, Lifshitz and Bonnor (for 
a review, see From there, one should expect that the conformal in- 

variance of YM lagrangian should result in a rather independent dynamics 
of the gauge fields regards to the expansion. The last should thus dilute the 
energy densities of the YM initial fluctuations. In fact, this is true when 
space-time is not significantly inhomogeneous i.e. when the perturbations 
of metric coefficients are negligible. In this case, the slightly inhomogeneous 
expansion mimics a conformal transformation of the metric. However, when 
the metric perturbations become important, we observe that the gauge fields 
start diffusing and the pressures and density contrasts start growing. 
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Before going further, let us examine more precisely the evolution of the gauge 
potentials. One can see in Figures El to |H1 that the gauge fields seem to roll 
on the potential well V = a 4 of the homogeneous solution but this potential 
is now distorded by the perturbations. For example, the electric component 
a, related to the speed of the oscillation, reaches its maxima at the middle 
of the phase space orbit a B = (see Figure EJ), while magnetic components b 
and c oscillate at almost the same frequency of the background solution o B ■ 
Let us give an idea of the distortion of the potential. Discarding the kinetic 
terms into Eq. (jl4j) . this potential can be written, in terms of the variables 
introduced in section 3.2: 



In Figure 01 we have represented the potential V associated to the results 
in Figures El 13 El dH and [T51 The oscillation has been unfolded along 
the a— axis and the resulting plot has been split in two for sake of readiness. 
The distortion has also been scaled up to be visible (the function plotted is 
a 4 + 6 x 10 3 (V — cr 4 ), the other parameters are same as in Figure EJ). The 
potential of inhomogeneous YM cosmologies appears to be an implicit func- 
tion of time. 

Furthermore, the gauge fields diffuse along the x~ ax i s > the perturbations 
spreading out with time. This diffusion is due to the increasing importance 
of spatial gradients of the gauge fields and the inhomogeneity of the metric 
in the YM equations. Another way to see this is to draw the evolution of 
the gauge fields during half a period of o B for the same amplitude of the 
initial perturbation but for different values of the background initial expan- 
sion rate (parametrized here by p B ). This is done in Figure ITUl for an initial 
perturbation in the electric component a. For small values of the expansion 
rate fFigure ITUl (i)). the diffusion is important but it rapidly disappears for 
stronger expansion (Figure ITUl (ii) and (in)). A much clearer way to express 
this is to compare the size of the fluctuations, given roughly by their coher- 
ence length Ap, to the Hubble length L H = l/H — t + J~-^ (c = 1). At the 
beginning of the evolution (t — 0), the last is given by 



Therefore, by varying the parameter p B and keeping the same coherence 
length through the same parameter w, we just modify the ratio between 




3 2/ 2 

V 



(36) 
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both quantities. The corresponding initial Hubble lengths have been spec- 
ified in the caption of Figure ITU1 with a coherence length Xp « 1. We can 
thus say that the diffusion is initially important for the short wavelengths 
Xp ~ Lh and, for wavelengths greater than the Hubble length, becomes im- 
portant only when the spatial inhomogeneity is strong enough, as in Figures 
li El El H and CU 

The observable quantities are represented in Figure ^2 under the form of 
the density contrast: 

P B 

and similar relations for the radial (T/) and tangential (T|) pressures con- 
trasts. It should also be noted that, according to the conformal invariance 
of YM fields, the scalar curvature TZ vanishes exactly and therefore the com- 
ponents of the stress-energy tensor are directly proportional to those of the 
Ricci tensor through Einstein equations and we have therefore a direct infor- 
mation on the spatial inhomogeneities. 

These contrasts oscillates in general, reflecting the underlying oscillations of 
gauge potentials, even though the amplitude of the oscillation is damped 
in the first stages of the evolution and then grow after some time, with a 
different regime for density and pressures contrasts. For example, in the 
case illustrated in Figure [TTJ the density constrast first decays, then enter a 
growing regime while the amplitude of the (radial and tangential) pressures 
contrasts slowly decreases. 

According to the radiation nature of YM, one should have expected that 
the evolution of the density contrast would have been a superposition of a 
growing mode in squared power of the conformal time 6 and a decaying one in 
inverse squared power of conformal time. This is well-known from a simple 
first order perturbation computation (see [25 ). In fact, when the perturba- 
tions of the geometry are small and their gradients negligible, the expansion 
can be considered as roughly homogeneous and it dilutes the energy of the 
conformally invariant YM fields. This case is a pretty good approximation 
of the early universe and is represented in Figure ^1 the amplitudes of the 
density and pressures contrasts decay with conformal time while the gauge 
potentials still oscillate in their distorded interaction potential, insensitive 
to an expansion that is too close to a conformal transformation. The grow- 
ing modes of the linear theory simply do not appear immediately due the 
scale invariance of YM fields. This means also that the YM cosmology can 

6 The density contrast is proportional to the world time in the synchronous gauge N = 1. 
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be much more inhomogeneous than it looks at the level of the observable 
quantities as the non observable ones (the metric components and the gauge 
potentials) do not decay in general, as we have seen before. 

By the way, when space inhomogeneity is too strong, growing modes for 
the contrasts and the perturbations of the gauge potentials (see Figure [T3] 
for illustration) appear, along with the diffusion already mentionned above. 
In fact, the details of this regime seem to be quite complex as it strongly de- 
pends on the interactions between the gauge components and on the field (s) 
initially perturbed. More work is under way to explore that regime and to 
characterize it in terms of the two dimensional abelian Higgs model of the 
theory. 

To summarize the evolution of Figures El El HI El HH if one starts with 
some fluctuations in a universe dominated by YM fields, the density and 
pressures contrasts will enter an oscillating regime, firstly damped by the 
quasi-homogeneous expansion. Then, when the expansion becomes signif- 
icantly inhomogeneous, the constrasts grow, at a time depending on the 
nature, the size and the amplitude of the initial perturbation. The transition 
to growing modes (and diffusion) seems indeed to occur when the perturba- 
tion of the metric coefficients 5 m = (m — 1)R or 5i = (/ — l)R 2 x 2 reach the 
order of magnitude of the scale factor (here a few percents, see Figure EJ). 
Therefore, the transition to a non conformal expansion regime is to be linked 
directly with the emergence of the instabilities. 

If the existence of growing modes is due here to the fact that the inho- 
mogeneous expansion looks no longer like a conformal transformation of the 
metric and therefore yields to non trivial interactions with the inhomoge- 
neous gravitational field, it would be interesting to see what happens to the 
growing modes with other gauges than the conformal one (N = R(t)) used 
here, for example, in the case of maximal slicing (cf. [21] and references 
therein). Indeed, this will lead to the study of a true gravitational collapse 
of YM fields instead of a shiny gravitational instability. 

Now, let us be more precise about how the parameters of the simulation 
affect the evolution of the observable quantities. A run is specified by the 
following parameters: pfi, the initial background density ; w, related to the 
coherence length X P of the fluctuations ; e, the relative amplitude of the 
perturbation ; er , the initial position of the Gal'tsov-Volkov particle in the 
interaction potential, and the nature of the perturbation (what fields were 
initially perturbed). The first two define the ratio between the coherence 
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and Hubble lengths jf-. Here, we mostly focus on perturbations larger than 
the Hubble length: in Figure ITTj X P > L H (Ap 3, L° H ps 0.09); in Figure 
C21 A P »> Lp: (A P ~ 2, L^- ~ 0.009) and finally in Figure CHI A P > L H 
(Ap ~ 2, L 1 ^ ps 0.3). In the first case, we see that the density contrast has 
grown after half a period of a ; in the second case and for the same period, we 
do not see yet any instability arising and finally, in Figure IT3~1 the contrasts 
have grown quite significantly during half a period of o. A more complete 
analysis shows that, all other parameters being equal, the large fluctuations 
Ap ^> L° H evolves more slowly than the small Ap ~ L° H . This can be intu- 
itively deduced from the previous comments on Figures ITT1IT3*1 and is clearer 
by examining Figure El where the evolutions of contrasts for two different 
coherence lengths are represented, all other parameters being equal. The evo- 
lution of the contrasts is qualitatively the same, although it is slower for large 
fluctuations (Figure ITU fa;), (v), (vi)). Please note the different duration and 
amplitudes of the two simulations. Therefore, the details of the evolution, 
e.g. the growing time, the period and/or the amplitudes of the contrasts, 
depend on the coherence length, the absolute amplitude of the perturbations 
(defined by J° for a certain field /) which is related to e and o"o, and the 

Jo 

nature of the perturbation. However, the scheme damped-then-growing oscil- 
lations seems to be general. We even performed some computations at short 
wavelengths Ap L H but we did not recover the usual Jeans criterion for 
radiation (cf. [25] )• We think that this criterion could be modified somehow 
due the locally anisotropic nature of YM pressures 7 (cf. [2E] and references 
therein for some properties of locally anisotropic fluids). More work is under 
way to fix that point. 

In Figure ^TJ (iv) , we see that the non-diagonal component of the stress- 
energy tensor Tq 1 becomes more and more homogeneous with expansion. This 
was observed in all simulations and can be deduced from Eqs. ([T7l - l33|) where 
one can see that this quantity evolves as the inverse fourth power of the scale 
factor R(t). 

For hyperbolic and closed cosmologies, we have performed some computa- 
tions that lead essentially to analogous features as those presented here in 
the flat case, for fluctuations of the order of the Hubble length. 

Finally, let us check out the violation of the three constraints TC, Tii and 

7 Locally anisotropic fluids are characterized by T| = Tg ^ in spherical symmetry 
and therefore have an equation of state different in the inhomogeneous case than the simple 
barotropic equation p = Kp with constant K . 
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Q s (see Figure The central bump in each constraint is related to the 
truncation errors when evaluating the spatial gradients, this can be improved 
by increasing the number of discretization points in the variable x- The in- 
teresting information is the time-evolution of the constraints: we see that the 
violation increases almost linearly during the computation and accelerates at 
the end. This is typical of the very simple explicit scheme we used here. 
However, it is important to note that the constraints remain several orders 
of magnitude below the perturbations of the related fields thanks to the very 
strong Courant condition we used. Another crucial point is the boundary 
effect as the careful reader will have already noticed from previous graphs. 
This can partly be improved by refining the accuracy of the IVP and partly 
by using more localized fluctuations. Some compactification of the \ domain 
could help in verifying more exactly the boundary conditions. In order to 
produce a more accurate method, one should also focus on improving the sta- 
bility by including implicit features in the scheme as well as a better control 
of boundary conditions. 

5 Conclusion 

By using a numerical method based on the hamiltonian formulation of the 
EYM system with SU(2) gauge group, we have checked that the primeval 
excitations of YM fields have been diluted by strong expansion in the very 
early universe, due to the scale invariance of YM equations. The gauge fields 
oscillate in a distorded potential well inherited from the homogeneous solu- 
tion while the metric suffers an ever increasing departure from homogeneity. 
Meanwhile, the corresponding density and pressures contrasts undergo firstly 
damped oscillations. 

However, the YM fields are gravitationally unstable in the long term. In- 
deed, the metric perturbations grow monotonically with time and, when the 
inhomogeneity of the metric and the gauge potentials cannot be neglected 
anymore, the last start diffusing non-linearly. This yields growing oscilla- 
tion modes for the density and pressures contrasts, with a transition from 
damping depending on the interactions between the perturbed gauge fields, 
the nature, size and amplitude of the primeval fluctuations. For example, 
the larger the fluctuation w.r.t. the Hubble length, the slower it will grow. 
A fruitful comparison can be made with the gravitational instability of the 
inflaton which lead us to the preliminary conclusion that the growth rate of 
the YM fluctuations should depend also on the strength of the gauge cou- 
pling constant and should be more important for strong interaction. 
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Nevertheless, in the real universe, we can conjecture that the phase tran- 
sitions that made YM fields short-ranged occur probably long before their 
inhomogeneities reach a sufficient size. Therefore such fields could be of mi- 
nor importance for cosmology. There are, however, two main objections to 
that affirmation. First, a more physical and quantitative analysis, including 
precise considerations about the amplitude and size of the primeval excita- 
tions of the gauge fields, other components of matter, precise cosmological 
parameters as well as quantum considerations on the status of the gauge 
fields after the phase transitions should be performed before concluding. In 
second, considering a Born-Infeld type modification of the YM lagrangian, 
deriving from string theory, which breaks the conformal invariance can be of 
interest. In particular, an interesting work would be to study the gravita- 
tional instability of non-abelian Born-Infeld cosmology which was examined 
in details in [TT] . 

Furthermore, another interesting point for cosmology is that the conformally 
invariant YM fields can hide for a while the deep inhomogeneity of space-time. 
Indeed, in a universe dominated by such fields, the observable quantities tend 
first to mimic a homogeneous universe while the metric is getting more and 
more inhomogeneous and the perturbations of the gauge fields oscillate. The 
question is to know which role such inhomogeneities could have played after 
the phase transitions that left the YM fields short-ranged as other types of 
matter fields feel differently the riddles of space-time hidden during the gauge 
dominating era than the conformally invariant gauge fields did. 

For all these reasons, we think that our results are of interest at least for 
the investigation of mathematical properties of the pure EYM theory, if not 
for a tiny step in the understanding of the very early universe . 
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Fi gUXe 1 '. Mean quadratic errors on the homogeneous solution for the gauge variables, (i) E a (ii) Ep 
(Hi) £ 7 (pg = 15, cr = 0, A x = 7 x 1CT 2 , At = 1Q- 7 ). 



SC.-G8 




Fig UX6 2i Mean quadratic errors on the homogi 
(same parameters as in Figure!?!). 




solution for the metric variables, (i) E m (ii) E\ 




Fi glire 3i Mean quadratic errors on the homogeneous solution for the constraints, (i) En (ii) E^i 
(Hi) Eg3 (same parameters as in Figure^. 
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Figure 4: Homogeneous solutions (i) R(t) (ii) a(t) (pg = 15, e c = 10" 4 , a = 0, w = 0.5, L° H a 0.09, 
Ax = 9.3 x 10" 2 , At = 1Q- 7 ). 




Figure 5i Perturbations of the metric coefficients, (i) (m — 1) R(t) (ii) (I — 1) R 2 x 2 (same parameters 
as in FigurelH) 





Figure Q'. Perturbation of the Erst gauge variable and the electric gauge potential, (i) (ct — <r) (ii) 
(a — Qfcfcg) (same parameters as in Figure^. 
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Fi gUXe 7'. Perturbation of the second gauge variable and the first magnetic gauge potential, (i) 
(fi — <r 3 ) (ii) (b — bi,)s 9 ) (same parameters as in FieurcTH. 




Fi glire 8: Perturbation of the third gauge variable and the second magnetic gauge potential, (i) 
(7 — <t 2 ) (ii) (c — Cjjfcg) (same parameters as in Figure^. 




Figure 9: Unfolded and rescaled shape of the potential V as a function of cr(t) and x ( see comments 
in the text). 
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Figure 10: Perturbations of the electric a (above) and magnetic b (middle) and c (below) potentials 
for different values of the initial expansion rate, (i) left: p^ = 0.15 (LP H f» 0.9) (ii) center: p^ = 1.5 
(L° H »J 0.3J (iii) right: pg = 15 (L° H « 0.09) (e a = 10~ 3 , w = 0.1, cr ( , = 0). 
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Figure 11: Density and pressures contrasts with non-diagonal component of the stress-energy tensor, 
(i) — jj- — 1 (ii) — — -g 1 (Hi) — '—g 1 (iv) Tg (same parameters as in Figure fill. 




Fig UX6 12i Example of decaying modes (early stages of the evolution), (i) (a — a B ) (ii) (b — b B ) (Hi) 
(c - c B ) (iv) (jfr - l) (v) (-2$. - lj (vi) - 1) (p B = 1.5 X 10 3 , L% » 0.009, e b = 10" 3 , 

w = 0.25, cr = 0, Ax = 0.15, At = 10~ 6 ). 
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